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O \ The aim of tliis study is to understand deeper the thermal diffusion transport 

process (Ludwig-Soret effect) at the microscopic level. For that purpose, the 
recently developed reverse nonequilibrium molecular dynamics method was 
used to calculate Soret coefficients of various systems in a systematic fashion. 
'^ . We studied binary Lennard-Jones (LJ) fluids near the triple point (of one of 

the components) in which we separately changed the ratio of one of the LJ 
O \ parameters mass, atomic diameter and interaction strength while keeping all 

other parameters fixed and identical. We observed that the magnitude of the 
Soret coefficient depends on all three ratios. Concerning its sign we found that 
p-^ . heavier species, smaller species and species with higher interaction strengths 

0^ , tend to accumulate in the cold region whereas the other ones (lighter, bigger or 

^^ \ weaker bound) migrate to the hot region of our simulation cell. Additionally, 

the superposition of the influence of the various parameters was investigated as 
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well as more realistic mixtures. We found that in the experimentally relevant 
parameter range the contributions are nearly additive and that the mass ratio 
often is the dominating factor. 
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I. INTRODUCTION 

Transport processes play an important role in the understanding of the prop- 
erties of liquid mixtures. Thermal diffusion (the Ludwig-Soret effect) is one 
of these intriguing processes. It characterizes the flux of matter in response 
to a gradient in temperature and has already been studied for almost 150 
years: The migration of atoms and molecules as a consequence of a temper- 
ature gradient was first reported by LudwigEl when studying sodium-sulfate 
solutions in 1856. Later, in 1879, Soreta observed the same effect in other 
electrolyte solutions. In the case of gases, it was predicted independently 
by EnskogS and by ChapmanQ, and later confirmed by the experiments of 
Chapman and DootsonEI. Further, several small organic molecule or poly- 
mer mixtures have been investigated by various experimental methods such 
as thermo-gravitational columns (e.g. Clusius and DickelEi), thermal field flow 
fractionation (e.g. GiddingsQ) or most recently by thermal diffusion forced 
Rayleigh scattering (KohlerB). There are numerous examples for the tech- 
nological significance of thermal diffusion, e.g. presented by Kincaid, Ratkje 
and IIafskj0ldH'li9. However, there is still no satisfactory theory to explain 
the effect. "It is the only hydrodynamic transport mechanism that lacks a 
simple physical explanation. "a. But the situation is even worse: One can- 
not even approximately predict the transport coefficients of closely related 
systems. The size of the Ludwig-Soret effect differs e.g. for equimolar mix- 
tures of benzene-chlorobenzene and benzene-nitrobenzene roughly by a factor 
of lOtil. Additionally, the experimental data base is rather small and partly 
inconsistentcl. One reason for this is probably that thermally driven flows of 
matter are minor effects and thus several orders of magnitude smaller than 
concentration driven flows. 

In multi-component systems, couplings between the different types of trans- 
port are, on a fundamental level, described by Onsager's linear relationshipsllj: 

N 

Jk = Y,^kjXj + LkgXq k = l,...,N (1.1) 

N 
Jq = 2_^ LqjXj + LqqXq 
J=l 

with Jfc being the flux of component k in gcm~-^s~^ with respect to the center 
of mass of the system and Jq being the heat flux in Jcm~ s~^. The Xa 
are thermodynamic forces, specifically for a binary mixture Xq = — ^^ and 
Xj = — VT(Aii — /U2)- (The subscript T denotes that the gradient has to be 
taken under isothermal conditions.) The L^p are the Onsager coefficients 
with a = {k, q) and N is number of species in the mixture. Note that this 
representation assumes J\\X, otherwise each L^/s has to be a 3 x 3 tensor. 
The assumption of linear response holds for many transport processes through 



condensed media. In this picture, thermal diffusion is characterized by the 
matrix coefHcient L^g. 



Although equations 1.1 have been known for several decades, there is yet no 
deeper insight into the microscopic mechanisms of the coupling of heat and 
mass transport. Liquid mixtures are often difficult to understand by analytical 
theories. Performing computer simulations, in contrast, is a way to look easily 
at the microscopic properties of a liquid system. The most promising results 
to understand thermal diffusion seem to come from molecular dynamics (MD) 
simulations. Several MD techniques have been developed in order to study 
the Ludwig-Soret effect. The early algorithms suffered from conceptual as 
well as practical problems: firstly, of how to define and compute the heat flow 
or the heat of transfer accurately in the microscopic picture and, secondly, 
from the large perturbation fields that were necessary to observe thermal 
diffusiontjilZI) . Hafskj0ld et al. E3 were the first to develop an algorithm 
that side-steps these problems. In fact, their methods are similar in spirit to 
ours. They define microscopic fluxes in an unambiguous way. However, they 
mostly investigated dilute gas mixturesu or interfacial systemsHScj so far. To 
our knowledge, there is no systematic study of thermal diffusion in the high 
density liquid regime up to present. 

In this contribution, we study both qualitative and quantitative aspects of 
the Ludwig-Soret effect. This is done by performing reverse non-equilibrium 
molecular dynamics (RNEMD) computer simulations for binary Lennard- 
Jones (LJ) systems. We investigate the influence of systematic variations of 
the physical parameters (mass m, atomic diameter a and interaction strength 
e) of our model system: By changing one parameter ratio (e.g. 7711/7712) while 
keeping all other parameters fixed and identical, we are able to observe ex- 
actly how the Ludwig-Soret effect depends on every one of them separately. 
We exploit the advantage of better control of the system variables in com- 
puter simulation over experiment to better understand thermal diffusion. We 
discuss our simulation results in the framework of linear irreversible thermody- 
namics. Moreover, we establish empirical rules for the influence of differences 
in the individual molecular parameters on the Ludwig-Soret effect and finally 
investigate to which extent they are additive. 



II. THEORY 

Linear Irreversible Thermodynamics 

In order to define proper quantities in the theory of linear irreversible ther- 
modynamics, the assumption of local equilibrium is essential. It enables us 
to apply the well-known equilibrium machinery to local volume elements of a 
perturbed system by defining intensive quantities derived from the extensive 
equilibrium quantitiesE^. The entropy source strength s given in —j^ — j is the 



fundamental quantity of this theory. Due to Onsager's variational principlellHI, 
s > has to be minimized. The non-equilibrium steady state is the one in 
which the least amount of total entropy is produced: 

s{f,t)dV — > min! (2.1) 

V 

It is also the state in which heat is most efficiently conducted through the 
system. Using a local form of the Gibbs equation in connection with the con- 
tinuity equation for the entropy (for details, see e.g. E3) leads to the following 
equation for s in a binary system with coupled flux of heat and matter: 

s{f, t) = - [4(r, t) ■ VT{f, t)] j;^ - [Ji • VrKr, t)] ^^ (2.2) 

with Ji being the mass flux of species 1 in —^ with respect to the center 
of mass of the system, Jq the heat flux in -^, fi = fii — fi2 the effective 
chemical potential in — and T the temperature in K. The right hand side 



of equation ^^ is the sum of products of the fluxes (of heat and matter), 
thermodynamic forces (the gradients) and of 'thermodynamic factors' J^- and 
7p, respectively. The latter provide for the correct dimensions. In this article, 
we consider the time-independent steady state in which the matter flux has 
died out and a constant temperature gradient is the only origin for the entropy 
source strength: 

^(-lj;=o = -[4(0-vr]^. (2.3) 

Combining this with Fourier's law for heat conduction we get: 

Kr1/,=o = -^- (^) • (2-4) 

The thermal conductivity depends on density, temperature and concentration 
of the mixture, which may vary over the system. Still, we assume that the 
perturbation and, hence, the variation in these quantities is small enough, so 
it is sensible to use an average thermal conductivity Xav The magnitude of 
the temperature gradient VT is the decisive quantity of the entropy source 
strength in our binary system. Note flnally, that the local entropy production 
increases towards the cold region of the system as 1/T(r)^. 



Phenomenological Transport Coefficients 



Equations 1.1 give the impression of a universal, clean and symmetric {La/3 = 
L^q) theory. Unfortunately, they have the disadvantage that the quantities 
of interest (i.e. the Onsager coefficients and the thermodynamic forces) are 



related to but not identical with experimentally measurable transport coef- 
ficients and fields. The relations between them and the Onsager coefficients 



Laj3 are well knownfcf. e.g. IHjc3). In their analysis, RNEMD simulations are 
akin to experiments, so in the following we will concentrate on the relevant 
experimental transport coefficients: the self diffusion coefficients D^ {k = 1, 2) 
(dimension err? j s) and the Soret coefficient St (dimension 1/i^). The former 
are calculated in our simulations via the Einstein route: 



6dt 



ff\t)-fj''\0)^). (2.5) 



Averaging is performed over time origins as well as over all particles j of type 
k. Since we empirically found that the D^ calculated out of RNEMD sim- 
ulations are (within statistical uncertainty) identical with equilibrium data 
(as we checked for some systems), we skipped further equilibrium runs to 
save resources and use the RNEMD data. The Soret coefficient is defined as 
St = Dt/Di2 and has the physical meaning of the relative strength between 
thermally induced diffusion (characterized by the thermal diffusion coefficient 
Dt) and concentration-driven, Fickian diffusion (characterized by the mu- 
tual diffusion coefficient D12). Phenomenologically, thermal diffusion is often 
expressed as 



Ji = -D12P 



dwi\ fdT 

+ Stwi{1 -wi) [ —- 



dz J \ dz 



Here, we are assuming field (temperature gradient) and fluxes (energy and 
matter) in z direction. Ji is the flux of species 1, p the average mass density 
(assuming that the temperature gradient and the resulting density gradient 
are small) and wi = ximi/{ximi + X2m2) the weight fraction of species 1 (x^ 



denotes the mole fraction of species k). Equation 2.6 appears to be similar to 
Onsager 's description (Equations ^]^). Note, however, that the gradient of the 
chemical potential has been replaced by the experimentally more accessible 
and technologically more relevant weight fraction gradient. If the system is 
continuously subjected to a temperature gradient, it will be driven to a non- 
equilibrium steady state: energy (heat) is then constantly flowing through 
it while the mass flux has stopped (Ji = 0) and a constant concentration 



In the special case of equimolar mixtures {xi = X2 = 0.5), equation 2.7 reduces 
further to 



«--(^)(S) 



(2.6) 



(or weight fraction) gradient has been established. Equation 2.(: can then be 
simplified and St, after application of the chain rule, be obtained from the 
temperature and molar fraction gradients in the system: 

9 = 1 f^\ f^y = 1 (dxi\ (dry 

^ wi{l-wi)\dz J \dzj xiil-xi)\dz J \dzj ■ ^'^ 



The gradients appearing in equation 2.8 can be calculated directly in our 



RNEMD simulations, making the analysis of the simulation data swift and 
easy. In particular, St can be calculated without previous evaluation of Dt 
and Di2- Throughout this work, positive values of St signify, that species 1 
tends to accumulate in the cold regions of the simulation box. 



III. COMPUTATIONAL DETAILS 

Model System and Simulation Details 

All simulations are performed with 1500 Lennard-Jones (LJ) atomic with 
cutoff distance re = 2.5(Ti, ui > a2- The potential is defined as: 



Ulj (r) = 4e 



(T\12 /(J 

r J \r 



We investigate equimolar mixtures of two species only. Consequently, our 
model system has six physical parameters: atomic masses nik {k = 1,2), 
atomic diameters (in form of Ufc) and the interaction strengths (in form of Sk) 
of our two model species. For the interaction between unlike particles, the 
Lorentz-Berthelot mixing rulesE3 are applied. No corrections to the interac- 
tion parameters due to excess functions (e.g. molar volume and enthalpy) are 
considered, although that can be done and might be useful in case of more re- 
alistic simulations. The orthorhombic periodic simulation cell is of size {L* x 
L* X 3-L*), with L* being roughly 8 — 13. The asterisk indicates Lennard-Jones 
reduced unitsEd, using L J- Argon values as reference values (Table 0). The vol- 
ume is held constant during the simulation. Our simulations are performed 
in the dense liquid state, i.e. p* = 0.8 - 0.85 and f* = 0.75 - 0.85 (Table |), 
which is close to the triple point of species 2. The average temperature T* 
is maintained at its value by the weak-coupling scheme of BerendsenE3 with 
a coupling time of t* = 4.66. We do so to prevent long-time drifts which 
might occur due to limited precision of the discrete, stepwise solution of the 
equations of motion. They are integrated using the velocity- Verlet algorithm 
in connection with a multiple-time-step schemec3 for the force calculations. 
Long-range forces (r > 1.7am with am = max{ai,a2}) are evaluated only 
every 4 time steps At* = 4.66 • 10^^, the switching range is 0.1am- Our sim- 
ulations are performed with run lengths of 3 • 10^ — 10 • 10^ steps. Shorter 
preparatory-runs of 0.8 • 10^ — 1.5 • 10^ steps are made to establish the non- 
equilibrium steady state. The runs have to be that long because the fluctua- 
tions in the mole fraction profile are relatively large and thermal diffusion is a 
weak effect compared to other transport properties. Positions and velocities 
are written out every 500 — 3000 steps, local temperature and concentration 
in the slabs (see below) every 50 — 100 steps for calculation of the relevant 
gradients. 
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(3.1) 



The RNEMD algorithm 

The RNEMD method reverses the usual cause-and-effect picture of non equi- 
librium simulations. For the study of thermal diffusion the "effect" , the heat 
flux, is imposed on the system whereas the "cause" , the temperature gradient 
is obtained from the simulation: We divide the simulation cell into A'^, = 20 
slabs perpendicular to the z-direction (Fig. ||). The slabs are chosen to be 
equally thick, i.e. to have identical volumes, and large enough so that reason- 
able statistics can be expected: each slab contains 75 particles on average. A 
flow of heat Jq is artificially maintained by exchanging velocities of suitably 
selected particles : Slab is defined as the 'cool' slab and slab A'^s/2 as the 
'hot' slab. Every N(,xch steps, we search through all atoms of the cold slab 
and determine the hottest ones of both species. In the hot slab, we proceed 
conversely and determine the coldest atoms of both species. Then we ex- 
change the velocities of the so determined atoms of the same species. We 
were always able to find at least one particle in the cold slab that is hot- 
ter than the coldest particle of the hot slab, since the Maxwell-Boltzmann 
distribution of atomic kinetic energies is very broad compared to the temper- 
ature difference of the two slabs. The exchange period has to be adjusted 
such that the perturbation is weak enough for linear response to hold (typ- 
ically: N^xch = 50 — 200 steps). This leads to a temperature variation of 
'^{Tlot - T:oM)/iTLt + T:,m) ^ 0.1 - 0.15. It was showJUli that for tem- 
perature differences of up to half the average temperature linear response 
is still approximately fulfilled. The unphysical velocity exchange leads to a 
physical flux of heat in the opposite direction through the intervening slabs. 
As the response is linear, the procedure causes linear profiles of temperature, 
overall density and concentration (Fig. ^). At steady state, physical and un- 
physical heat fiux have the same magnitude because of energy conservation. 



Thus, cf. equation 2A, the temperature gradient will be minimal. To avoid 
a possible breakdown of the local equilibrium due to the unphysical energy 
transfer in the thermostatting slabs, these are excluded from gradient cal- 
culations. General properties and details of the RNEMD technique as well 
as how to cornnute other transport coefficients with it have been published 
elsewhereE3'E3'E3. 



IV. RESULTS AND DISCUSSION 

Data Comparison with other Work 

Since the Ar-Kr LJ-system was mostly investigated by MD computer sim- 
ulations in the past, we made a consistency check of our results with such 
a system. Table n compares our own results with data obtained by several 
other groups, all using identical LJ-parameters for Ar and Kr, respectively 



(Table pl[ ). For the self diffusion coefficients, the agreement is excellent. The 
situation is different concerning the Soret coefficient. It can be seen, that 
although the signs for St are identical and all values lie close to each other, 
they do not agree with each other within the statistical uncertainty. However, 
comparison with the older results is somewhat difficult since system sizes were 
much smaller (128 or 256 particles). Moreover, we found in our simulations 
that the value of St is strongly influenced by the specific state point at which 
the simulation was run. Minor shifts in the density were e.g. sufficient to 
change the value of St by more than 30% (Table O) . Vogelsang and Hoheisel 
reportedtj, that their values are systematically too low due to an inaccurate 
determination of the partial enthalpies. Extrapolation to zero field strength 
(MacGowan and EvansuO, Paolini and CicottitJ) or zero temperature gra- 
dient (as in our case) are also sources of error. Figure ^ shows that in our 
simulations the Soret coefficients can be well extrapolated to zero field. The 
two outermost points deviate from the extrapolation line for two different 
reasons. At low temperature gradients (left) there is too much statistical un- 
certainty in the concentration gradient. At high gradient (right) the linear 
response breaks down because of the too large perturbation. Additionally, 
the temperature of the cool region lies then below the freezing temperature. 



Systematic Variations of Parameters 

As long as we change only one parameter ratio at a time, the problem is 
symmetric in species identity, i.e. it is sufficient to consider ratios larger than 
1.0. We define species 2 as reference species and identify it with L J- Argon: 
7712 = 39.95 amu, cr2 = 0.3405 nm and £2 = 1-0 kJ/mol. When changing 
the mass ratio while keeping all other parameters fixed and identical, the 
positive sign of St indicates that the heavier species 1 prefers the cold side of 
the simulation box whereas the lighter species 2 favors the hot side (Fig. |^. 
Moreover, St rises monotonically for the whole mass ratio range. As a side 
result we found that the behaviour of the self diffusion coefficients could be 
very well fitted by power laws: Di oc (7771/7/1,2 )~^'^^^ and D2 oc (777i/77i2)~*''^^^, 
which is very close to exponents —3 and —j. They are also similar to those 
first found by Bearman and Jolly in simulations of Ar-Kr mixtures (—0.3 
and —0.2, respectively)^^. For the interaction parameter e it can be seen 
that St is positive for £1 > £2 (Fig. ^), i.e. species with deeper potential 
wells prefer the cold side of the simulation box while the more weakly bound 
species accumulate on the hot side. The effect of mass and interaction strength 
variations is intuitively understandable since both times, the species with the 
lower mobility (i.e. lower self diffusion coefficient, cf. Fig. ^, |5|) favors the cold 

For the atomic diameter parameter a, in contrast, we find that St is nega- 
tive for o"! > (T2, i.e. smaller particles prefer the cold side of the simulation 
box although their mobility is higher compared to the bigger species (Fig. |^. 
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Furthermore, a monotonic behavior was only found for ai/(T2 < 1.25 (which 
would already be large for realistic mixtures). A pressure effect can here be 
ruled out as the pressure did not vary more than 5%. For higher ratios, the St 
data points do not show a monotonic behaviour any more. Here, we also found 
a large disparity of the two self diffusion coefficients. Specifically, the bigger 
species becomes increasingly immobile compared to the smaller species and 
the absolute number of bigger particles per slab is nearly constant (Fig. IW). 
That indicates that a regime is reached in which the mobility of the larger 
particles is massively hindered by their own space-filling arrangement in the 
box. Therefore, the large particles can no longer exhibit a density gradient. 
As a consequence, the small species remains the only one capable to accumu- 
late in the cold region: It is done for entropic reasons (the cooler a region, the 
higher the overall density) . Now the question arises if the trend 'the bigger the 
species, the more it tends to accumulate in the hot region' still holds if both 
species can move around unhindered. Therefore, a run with lower density and 
higher temperature (o"i/o"2 = 1.4, p* = 0.6, T* = 1.15) was performed. The 
bigger species then builds up a concentration gradient and accumulates at 
the hot side with 5"^ = — 6.1 ±0.6, confirming the trend observed for atomic 
diameter ratios. However, the high ratio regime is artificial and the above ex- 
planation not applicable in the realistic area of (Ti/(T2 < 1.25. There, another 
argument seems to be more appropriate: If we combine the LJ parameters a 
and e into the variable e = -% we obtain a measure for the stored energy per 
volume, a potential or cohesive energy density. The Soret coefficient varies 
linearly with the ratio ei/e2 (Fig. ^). The species with higher cohesive energy 
density accumulates in the cold region of the system. 



Realistic systems and superposition of different contributions to the Soret coefficient 

We also investigated how more complicated mixtures (species differing in all 
parameters) behave under the influence of a temperature gradient. Thus, we 
can check if the parameters contribute more or less independently to the Soret 
coefficient and if one of them is dominating. The investigated systems com- 
prised all possible liquid binary mixtures of Ar,Kr,Xe and CH4 (Table ^, 
representing species with strongly deviating parameters. The following em- 
pirical laws were obtained by low order fits of the independent parameter 
variations (Table ^Vj): 



r'/K] = -0.7-(^y + 9.5.(^ 



SrilO-yK] = -0.7 ■ [^] + 9.5 • ^ ^ - 8.8 for mi/ms ^ 8.0 (4.1) 
\m2j \m2j 

St[W-^/K] = 67.4 • ( — J - 179.3 • [ — J + 111.9 for cJi/cJa < 1.25 (4.2) 

St[IQ-^/K] = 4.4 • f — j + 3.5 • f — j - 7.9 for ei/ej ^ 1.75 (4.3) 



Using the additivity of the three contributions, we predicted the Soret coeffi- 



cients of the reahstic systems and compared them, using Eqns. |4 .11 , 14.21 and |4.3| , 
with the simulated values (Fig. |9|). The agreement is excellent: The sign of St 
is correct for all cases and the relative deviation is below 30% for all realistic 
mixtures. The results show also clearly that increasing Soret coefficients are 
mainly correlated with increasing mass ratios (Fig. |9|). That is because in 
realistic systems (of unlike noble gases), the mass difference between the two 
species is usually much larger than the deviations in their effective diameter 
or interaction strength. Additionally, the effects of a and e tend to cancel 
out each other, since mostly, the bigger atom is also the one with the deeper 
potential well. For this reason it is possible to combine the last two equations 
(|4.2| and |4.3|) into a single one which involves the cohesive energy densities. 



That is done by fitting a straight line to the composed data shown in Fig. ^: 

St[W-^/K] = 14.8 • ( — J - 14.8 (4.4) 



As it can be seen in Table ^, the accuracy of the prediction using Eqns. 4.1 



and 4.4 is the same compared to Eqns. 4.1, ^^ and 4.3 for most data points. 
Only in the case of the Xe-Ar mixture, the predicted values differ significantly. 
Here, the latter prediction is much lower than the first one. That is most 
probable because the individual parameter ratios are large (and, hence, tend 
to be more inaccurate) while, due to compensation effects, the composed 
ratio does not differ much from unity compared to the other mixtures. To 
check if the empirical laws are applicable to experimental data, we compared 
the predicted results for equimolar mixtures of benzene-chlorobenzene and 
benzene-nitrobenzene (mentioned in the introduction) to the experimental 
values. The sign as well as the order of magnitude of St could be reproduced: 
For benzene-chlorobenzene, S'T(pi'ed) = —5.0-10"^ (exp.: —1.1 -10"^) and for 
benzene-nitrobenzene /S'T(pred) = —12.6 • 10~^ (exp.: —11.0 • 10~^). In view 
of the simplicity of the model, we consider this as an excellent result which 
motivates us to continue work in this direction. 



V. CONCLUSIONS 

We investigated with computer simulations thermal diffusion processes in sim- 
ple two-component Lennard-Jones model liquids out of equilibrium. We ob- 
served that for dense liquid mixtures, the size of the Ludwig-Soret effect, 
depends on the ratios of all three model parameters: the mass, the atomic 
diameter and the interaction strength. However, the combination of the two 
latter into a single parameter e = -%, the potential energy density of a species, 
seems to be especially appropriate to explain the sign of the Ludwig-Soret ef- 
fect. 

Parts of our results are in line with earlier invest iga tions, especially by Vo- 
gelsang and Hoheisel et aLE^'ES and Hafskj0ld et al.E3a, while others are com- 

10 



pletely new. For isotopic systems, the lighter species can transport a higher 
amount of heat through the system if it accumulates in the hot region. More- 
over, the possibility to superpose the influence of the various parameters was 
investigated on more realistic mixtures of LJ noble gases. For most of these 
systems, the mass ratio is the dominating factor that determines sign and 
value of St- This trend is for realistic mixtures augmented by the fact, that 
the effects of a and e tend to cancel out each other, since bigger species often 
also are the ones with the deeper potential depth. Further investigations will 
compare our simulation results with experimental data. 
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TABLES 



Quantity Reduced units 

Mass m* = ^ 



Length 



r 

CFAr 



E 



Energy E* 



Time i* = ^ J^^ 



Diffusion coefficient D* = -D_ mAn 

<^Ar y SAr 

Soret coefficient S^ = Stj^^ 

Pressure p* = p-^ 

^Ar 



Density p* 

Temperature T* 



Ml 
V 



kgT 
£12 



TABLE L Reduced units for physical quantities. LJ-Argon was taken to be the reference 
system (Table pl| ). The mixed quantities ai2 and £12 are derived by the Lorentz-Berthelot rules. 
For density and temperature they were used to conserve the state point in the phase diagram of 
the mixed system. 
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Ref. 


Method^ 


Number 
of atoms 


P* 


T* 


St 
[10-3^-1] 


DAr 

[lO-^cm^s-i] 


DKr 

[lO-^cm^s-i] 


ME2 


NEMD 


108,256 


0.7902 


0.805 


11.3 






PC^ 


NEMD 


108,256 


0.7902 


0.824 


16.2±2.0 






VH^ 


EMD 


108,256 


0.803 
0.79 


0.81 
0.80 


9.1±2.0 


2.97±0.08 


2.60±0.08 


SRS 


EMD 


256 


0.79 


0.81 




2.97±0.08 


2.44±0.08 


this work 


NEMD 


1500 


0.797 


0.805 


10.5±1.3 


2.98±0.05 


2.48±0.05 


this work 


NEMD 


1500 


0.81 


0.85 


14.4±1.2 


2.97±0.05 


2.47±0.05 



TABLE II. Comparison of calculated values for the Soret and self diffusion coefficients of a 
Lennard-Jones Ar-Kr binary mixture. 

^ EMD: Equilibrium Molecular Dynamics, NEMD: Non-Equilibrium Molecu- 
lar Dynamics 
^ MacGowan and EvansEj 
3 Paohni and CicottiEJ 
^ Vogelsang and HoheisellB 
^ Schon and HoheiselcJ 
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Atom species 



m\amu\ 



m 



a\nm\ 



a 



e[kJ/mol] 



e^ 



Ar 

Kr 

Xe 

CH4 



39.95=1.00 

83.80=2.10 

131.29=3.29 

16.04=0.40 



0.3405=1.00 
0.3633=1.07 
0.3975=1.17 
0.3740=1.10 



1.00=1.00 
1.39=1.39 

1.72=1.72 
1.27=1.27 



TABLE III. Lennard-Jones parameters for the atoms and molecules of this study. Molecules 
are treated as single Lennard-Jones atoms. 
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rai/ra2 


0-1 /0-2 


£1/^2 


St 
[10-3^-1] 


[lO-^cm^s-^] 


D2 

[lO-^cm^s-^] 


P 

[lO^Pa] 


1.2 


1.0 


1.0 


1.4± 0.2 


2.96 


2.96 


4.6 ± 0.4 


1.5 


1.0 


1.0 


5.0± 0.5 


2.79 


2.82 


4.6 ± 0.3 


2.0 


1.0 


1.0 


8.9± 0.7 


2.63 


2.69 


4.7 ± 0.3 


3.0 


1.0 


1.0 


15.7± 0.5 


2.30 


2.37 


4.7 ± 0.3 


4.0 


1.0 


1.0 


19.4± 0.9 


2.09 


2.31 


4.6 ± 0.3 


6.0 


1.0 


1.0 


23.4± 0.6 


1.69 


1.96 


4.6 ± 0.4 


8.0 


1.0 


1.0 


24.6± 0.8 


1.51 


1.79 


4.6 ± 0.4 


1.0 


1.0 


1.1 


1.2± 0.3 


2.99 


3.00 


4.9 ± 0.4 


1.0 


1.0 


1.2 


2.7± 0.4 


3.03 


3.06 


5.0 ± 0.4 


1.0 


1.0 


1.3 


4.1± 0.5 


3.10 


3.23 


5.2 ± 0.4 


1.0 


1.0 


1.4 


5.5± 0.8 


3.12 


3.34 


5.3 ± 0.4 


1.0 


1.0 


1.5 


7.3± 0.7 


3.14 


3.40 


5.3 ± 0.4 


1.0 


1.0 


1.75 


11.6± 0.6 


3.17 


3.53 


5.5 ± 0.4 


1.0 


1.0 


2.0 


21. 0± 0.4 


3.15 


3.71 


5.8 ± 0.5 


1.0 


1.0 


2.25 


30.7± 0.9 


3.13 


3.87 


5.9 ± 0.6 


1.0 


1.0 


2.5 


41.4± 2.6 


3.05 


4.00 


6.0 ± 0.7 


1.0 


1.05 


1.0 


-2.6± 0.5 


2.90 


3.11 


4.1 ± 0.3 


1.0 


1.1 


1.0 


-3.6± 0.5 


2.89 


3.18 


3.9 ± 0.3 


1.0 


1.15 


1.0 


-5.4± 0.5 


2.78 


3.21 


3.9 ± 0.3 


1.0 


1.2 


1.0 


-6.1± 0.5 


2.62 


3.27 


3.8 ± 0.2 


1.0 


1.25 


1.0 


-7.1± 0.5 


2.56 


3.31 


3.8 ± 0.2 


1.0 


1.3 


1.0 


-6.1± 0.6 


2.46 


3.29 


4.0 ± 0.2 


1.0 


1.35 


1.0 


-5.7± 0.6 


2.37 


3.17 


4.2 ± 0.2 


1.0 


1.4 


1.0 


-7.0± 0.5 


2.09 


3.16 


4.4 ± 0.2 


1.0 


1.5 


1.0 


-5.2± 0.5 


1.72 


3.04 


5.0 ± 0.2 


1.0 


1.7 


1.0 


-8.4± 0.7 


1.14 


2.45 


6.2 ± 0.1 


1.0 


1.9 


1.0 


-9.6± 1.7 


0.58 


1.84 


7.6 ± 0.1 



TABLE IV. Thermal Diffusion data: Separate variation of every species 1 parameter at the 
state point T* = 0.85, p* = 0.81. Species 2 corresponds to Lennard- Jones Argon. The error is 
estimated by hnear regression (in case of St of the temperature and mole fraction profiles). For 
the self diffusion coefficients, the error is generally ±2%. 
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Model 

system 


mi/m2 


<7l/o"2 


£1/^2 


ei/e2 


simulated 

St 
[10-3^-1] 


predicted^ 

St 
[10-3^-1] 


predicted^ 

St 
[W-'^K-^] 


absolute^ 
deviation 
[lO-^K-^] 


relative^ 
deviation 

[%] 


An,, - Ar 


1.50 


1.00 


1.00 


1.00 


5.0±0.5 


3.9 


3.9 


1.1 


22.0 


Xe-Kr 


1.57 


1.03 


1.23 


1.13 


5.1±0.7 


6.2 


6.3 


-1.1 


-21.6 


Kr - Ar 


2.1 


1.07 


1.39 


1.13 


14.4±1.2 


10.8 


10.1 


3.6 


25.0 


Ar - CHi 


2.49 


0.91 


0.79 


1.05 


9.3±0.7 


10.7 


11.2 


-1.4 


-15.0 


Xe-Ar 


3.29 


1.17 


1.72 


1.07 


18.6±0.8 


20.4 


16.0 


-1.8 


-9.7 


Kr - CHi 


5.22 


0.97 


1.10 


1.21 


22.3±0.6 


24.3 


24.8 


-2.0 


-9.0 


Xe - CHi 


8.17 


1.06 


1.35 


1.13 


23.1±0.7 


24.5 


24.1 


-1.4 


-6.1 


Ar^,, - Ar 


8.00 


1.00 


1.00 


1.00 


24.6±0.8 


22.4 


22.4 


2.2 


8.9 



TABLE V. Soret coefficients for several simple atomic systems at the state point T* = 0.85, 
p* = 0.81. The run length was 4.0 • 10^ steps. Additionally, the data for two isotopic systems 
is presented. A comparison indicates, that St is chiefly governed by the mass ratio. The error 
is estimated by error propagation from the errors of the slopes of temperature and mole fraction 
profiles as obtained in the linear regression. 



^ using Eqns. 4.1,4.2 and 4.3 
^ using Eqns. 4.1 and 4.4 
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FIG. 1. Subdividing the periodic simulation box into slabs. Slab is defined to be the cool 
slab, slab Ns/2 to be the hot slab. Kinetic energy is artificially transferred by the heat exchange 
algorithm and flows back by thermal conduction. 






X) 




c^ 


78 


K/^ 




^ 


76 






C/D 


74 


<U 




O 




• ^H 


72 


■t— 1 




^ 




Dh 



0123456789 10 
slab number 

FIG. 2. Typical profiles of system properties at the state point T* = 0.81, p* = 0.85. 
ci/(72 = 1.2 and species 2 corresponds to Lennard-Jones Argon. Slab is defined to be the 
hot slab, slab 10 to be the cool slab. 
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AT in K / slab 



FIG. 3. Linear response regime and extrapolation to zero temperature gradient. Model system: 
mi/m2 = 4.0, o"i/iT2 = 1.0 and £1/62 = 1-3 where species 2 corresponds to Lennard-Jones Argon. 
A straight line was fitted to the data points (with the highest and the lowest point omitted) in a 
linear scale. The extrapolation value is St = 22.3 • 10~^/i^. 
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FIG. 8. The Soret coefficient varies approximately linearly with the ratio of 'cohesive energy 
densities' ei/e2 > l(efc = efc/o"|). Circles: The data for ei/e2 < 1 were calculated from 0"i/cj2 
variations at constant ExJEi = 1. For the data at ei/e2 > 1, ei/e2 was varied at constant oxja^ = 1- 
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